
library (foreign)
library (survival)
library (arm)
library (pscl)
library(apsrtable)


######################################################
#
#   Negative Binomial for MID Initiation
#   Thursday, July 26 1:55 p.m.
#######################################################




## Laptop   

TermLimitsData<-  read.dta("/Users/Jeff/Dropbox/TermLimits/Replication/Data/TermLimitsAnalysisData.dta")

summary(TermLimitsData)
#TermLimitsData <- subset(TermLimitsData, year >= 1960)
#attach(TermLimitsData)


##  Model   ##

mod1<- glm.nb (countinit ~ lameduck + hawkish + 
                           rivalry + numbord + parliament + gender +
                                   peace + peace2 + peace3,data=TermLimitsData)
summary(mod1)
odTest(mod1)
logLik(mod1)



#> odTest(mod1)
#Likelihood ratio test of H0: Poisson, as restricted NB model:
#n.b., the distribution of the test-statistic under H0 is non-standard
#e.g., see help(odTest) for details/references

#Critical value of test statistic at the alpha= 0.05 level: 2.7055 
#Chi-Square Test Statistic =  3.5144 p-value = 0.03042 
#> logLik(mod1)
#'log Lik.' -807.3266 (df=11)




forms <- strsplit(as.character(mod1$call)[2], split="~", fixed=T)
forms <- unlist(lapply(forms, function(x)paste("~", x, sep="")))
count.form <- as.formula(forms[2])

pred.dat <- data.frame(
  lameduck = c(0,1), 
  hawkish = c(0,0),
  rivalry = c(0,0),
  numbord = c(2,2),
  parliament = c(0,0),
  gender = c(0,0),
  peace = c(5,5),
  peace2 = c(25,25),
  peace3 = c(125,125)
) 




X.count <- model.matrix(count.form, data=pred.dat)

#all.b <- do.call("c", mod1$coef)
set.seed(22)
b.mat <- mvrnorm(1000, mod1$coef, vcov(mod1))

count.preds <- exp(X.count %*% t(b.mat))

count.ci <- t(apply(count.preds, 1, quantile, c(.5,.025,.975)))





pred_civ_account<-mean(count.preds[1,])
pred_civ_account_low <- count.ci[1,2]
pred_civ_account_high <- count.ci[1,3]


pred_civ_tl<-mean(count.preds[2,])
pred_civ_tl_low <- count.ci[2,2]
pred_civ_tl_high <- count.ci[2,3]





diff_civ <- count.preds[2,]-count.preds[1,]


diff_civ_mean <- mean(diff_civ)
diff_civ_low <- quantile(diff_civ, 0.025)
diff_civ_high <- quantile(diff_civ, 0.975)




#########################################
## Generating Graphing Data Set
#########################################


PredictedValues <- data.frame( pred_civ_account,pred_civ_account_low,pred_civ_account_high,
                                                                         pred_civ_tl,pred_civ_tl_low,pred_civ_tl_high,
									diff_civ_mean,diff_civ_low,diff_civ_high)

write.dta(PredictedValues,file="/Users/Jeff/Dropbox/TermLimits/Replication/Models/IndexEmpirical/Unconditional/MultipleDisputes/PredictedValues.dta")



